close all
clear
clc

% problem: Ur=2 and Er=6, 1ft(30.48cm) think slab, calculate reflectance and
% transmittance from 0~1Ghz

% Common value for simulation
c0 = 299792458; % light speed
gigahertz = 1e9;
ER_slab = 6;
UR_slab = 2;
ER_src = 1; % ER at source
UR_src = 1; % UR at source
nmax = sqrt(ER_slab*UR_slab); % refractive index for material
nbc = 1; % refractive index at the grid boundaries
nsrc = 1; % refractive index at the soruce cell
freq_max = 1*gigahertz; % 1GHz
slab_length = 30.48/100; % 1ft=30.48cm

% =================================================================
% Start calculate grid resolution
NRES = 20; % resolve the wave with at least 10 cells, better >=10
lambda_min = c0./freq_max./nmax;
delta_lambda = lambda_min/NRES;

NDRES = 4; % normally 1~4 for resolution of feature size
delta_d = slab_length/NDRES;

dz_temp = min(delta_lambda, delta_d);
Nz_slab = ceil(slab_length/dz_temp); % get the cell number for slab
dz = slab_length/Nz_slab; % get the size of grid

space_size = 10; % 10 cell at each side as space at least
% Calculte total Nz = slab:Nz_slab; + space=10 cell/each side; + 3 cells
% extra 3 cells:  1 for record reflection(Nz=1); 1 for TF/SF source(Nz=2);
%                 1 for record transmission(Nz=end, last cell)
Nz = Nz_slab + 2*(space_size) + 3;
% End calculate grid resolution
% =================================================================

% Initialize Materials to free space
ER = ones(1,Nz);
UR = ones(1,Nz);

% starting cell of slab
nz1 = 2 + space_size +1; % 1cell reflection+1cell sf/tf source + space_size
% end cell of slab
nz2 = nz1 +  Nz_slab -1; % nz1 + cells of slab -1(last cell of slab)

ER(nz1:nz2) = ER_slab;
UR(nz1:nz2) = UR_slab;

% =================================================================
% Start calculate delta_t and tau 
dt = nbc*dz/(2*c0);

tau = 1/(2*freq_max);
t0 = 6*tau;
t_prop = nmax*Nz*dz/c0;

T_total = 12*tau + 5*t_prop;
STEPS = ceil(T_total/dt);
% =================================================================

% =================================================================
% Start calculate Source for Ey/Hx mode
t=(0:STEPS-1).*dt; % time axis
delta_t = nsrc*dz/(2*c0) + dt/2; % total delay between E and H
A = -sqrt(ER_src/UR_src); % amplitude of H field
Esrc = exp(-((t-t0)/tau).^2); % E filed source
Hsrc = A*exp(-((t-t0+delta_t)/tau).^2); % H field source
% =================================================================

% =================================================================
% Inilize fourier transforms
NFREQ = 100;
freq_fft = linspace(0,freq_max,NFREQ);
K=exp(-1i*2*pi*dt.*freq_fft);
REF_fft = zeros(1,NFREQ); % reflection
TRN_fft = zeros(1,NFREQ); % transmission
SRC_fft = zeros(1,NFREQ); % srouce

% record time domain for compare
REF_t = zeros(1,STEPS);
TRN_t = zeros(1,STEPS);
SRC_t = zeros(1,STEPS);
% =================================================================

% Compute updated coefficients
mEy = (c0*dt)./ER/dz;
mHx = (c0*dt)./UR/dz;

% Initialize Ey and Hx to zero
Ey = zeros(1,Nz); 
Hx = zeros(1,Nz);

E3=0; E2=0; E1=0;
H3=0; H2=0; H1=0;

fig=figure;
set(fig,'Name', 'FDTD 1D SF/TF Simulation');
set(fig,'NumberTitle', 'off');

% position of the source, 2nd cell
Nz_src = 2; % 1st cell is for reflection record, 2nd cell for inject source
% Main FDTD Loop
for T = 1 : STEPS
    
    % Update H from E
   for nz = 1 : Nz-1
       Hx(nz) = Hx(nz) + mHx(nz)*(Ey(nz+1) - Ey(nz));
   end
   % Absorting Boundary Conditions
   Hx(Nz) = Hx(Nz) + mHx(Nz)*(E3 - Ey(Nz)); 
   
   % Handle H source, update SF/TF in H field
   Hx(Nz_src-1)=Hx(Nz_src-1)-mHx(Nz_src-1)*Esrc(T);
   % handle H at boundary
   H3=H2; H2=H1; H1=Hx(1);
   
   % update E from H
   % Absorting Boundary Conditions
   Ey(1) = Ey(1) + mEy(1)*(Hx(1) - H3); 
   for nz = 2 : Nz
       Ey(nz) = Ey(nz) + mEy(nz)*(Hx(nz) - Hx(nz-1));
   end
   % Handle E source, update SF/TF in E field
   Ey(Nz_src) = Ey(Nz_src)-mEy(Nz_src)*Hsrc(T);
   % handle E at boundary
   E3=E2; E2=E1; E1=Ey(Nz);

   % update fourier transforms
   for nf = 1:NFREQ
       REF_fft(nf) = REF_fft(nf) + (K(nf)^T)*Ey(1); % check reflection in k=1
       TRN_fft(nf) = TRN_fft(nf) + (K(nf)^T)*Ey(Nz); % check transmission in k=Nz
       SRC_fft(nf) = SRC_fft(nf) + (K(nf)^T)*Esrc(T); 
   end

   REF_t(T) = Ey(1);
   TRN_t(T) = Ey(Nz);
   SRC_t(T) = Esrc(T);
   
   %Plot
   imagesc(1:Nz,0,ER);
   % colormap('jet');
   hold on
   
   plot(Ey,'-b','LineWidth',2);
   hold on
   plot(Hx,'-r','LineWidth',2);

   title(sprintf('Step: %d of %d',T, STEPS));
   xlim([1 Nz])
   ylim([-1.5 1.5])
   drawnow;
   hold off
     
end

% move this out of loop to speed up for fourier transform calculate
REF_fft = REF_fft * dt;
TRN_fft = TRN_fft * dt;
SRC_fft = SRC_fft * dt;

TRN = abs(TRN_fft./SRC_fft).^2;
REF = abs(REF_fft./SRC_fft).^2;
CON = REF+TRN;
plot(freq_fft, TRN)
hold on
plot(freq_fft, REF)
hold on
plot(freq_fft, CON)








